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GEOPHYSICS,  should  the  paper  be  accepted  for  publication. 


A procedure  is  developed  for  obtaining  max iroun- likelihood  estimates 
of  the  reflection  coefficient  sequence  from  seismic  data.  The 
reflection  coefficient  sequence  is  modeled  as  an  impulsive  process, 
where  the  reflection  locations  are  statistically  independent  and  the 
reflection  amplitudes  are  uncorrelated  and  Gaussian  distributed.  The 
wavelet  and  all  other  parameters  are  assented  known.  The  results  of  this 
procedure  are  danonstrated  for  synthetic  data. 


I.  introduction 

The  optimal  smoother  method  for  estimating  reflection  coefficients 
from  seismic  reflection  signals,  as  developed  by  Mendel  [1-4],  is  a 
minimus-variance  linear  estimator  which  does  not  take  into  account  the 


impulsive  nature  of  the  reflection  coefficient  sequence.  Figure  1 shows 
the  results  of  using  the  optimal  another  on  a synthetic  seianogram;  the 
circles  mark  the  non- zero  reflection  coefficients  and  the  bars  depict 
the  corresponding  estimates.  As  one  can  see,  the  minimus-variance 
estimates  are  generally  non-zero  at  every  time  point  and  tend  to  under- 


□ □ 


As  with  the  minimum- variance  estimator,  we  model  the  seismic  signal 
as  a linear  convolution  of  the  source  wavelet  with  a reflection  process, 
corrupted  by  various  recording  effects.  These  recording  effects  include 
the  geophone  response,  aliasing  filter,  instrunent  noise,  cultural 
noise,  and  can  include  modeling  errors.  Further,  the  signal  to  be 


the  sane  as  the  signal  to  be  convolved  with  the  wavelet,  namely  the 
impulse  response  of  the  Barth.  The  difference  between  these  two  signals 
can  be  regarded  as  geological  effects,  which  include  spherical 
divergence,  absorbtion,  and  multiples.  Sane  of  these  effects  can  be 
included  in  our  model,  which  is  shown  in  Figure  2.  (Including  multiples 
would  require  a non-linear  model  of  very  high  order,  and  is  generally 
not  feasible  at  this  point.)  Given  this  model  in  state  vector  form,  and 
the  variance  of  y(k) , one  can  perform  optimal  smoothing. 


series  of  distinct  reflectors,  so  that  the  reflection  coefficient 


sequence  will  consist  of  a few  isolated  impulses.  To  avoid  a lengthy 
discussion  for  conversion  from  continuous  to  discrete  time,  let  us 
assune  that  these  reflections  can  only  arrive  at  discrete  times.  Due  to 
the  difficulty  in  modeling  multiple  reflections  and  transmission 
effects,  let  us  simply  ignore  the  additional  information  available  from 
these  effects.  Finally,  we  assune  that  the  locations  (in  time)  of  the 
reflections  are  statistically  independent  of  each  other,  and  that  the 
magnitudes  of  these  reflections  are  uncorrelatad  and  Gaussian 


Figure  3 depicts  a synthetic  reflection  coefficient  sequence 


generated  using  this  model,  and  Figure  4 depicts  the  corresponding 
seianograo.  For  this  example  we  have  ignored  all  geological  effects  and 
all  recording  effects  except  for  correction  by  additive  white  Gaussian 
noise.  This  is  not  a constraint  on  the  method,  but  was  done  simply  to 
avoid  complications  in  the  presentation.  Figure  4 depicts  the  signal 
from  which  the  estimates  in  Figure  1 were  obtained. 


II.  The  Product  Model 

To  simplify  the  estimation  problem,  we  have  found  it  convenient  to 
model  the  reflection  coefficient  sequence,  y(k),  as  a product  of  two 
random  processes  of  the  form 

w(k)  - r(k)  q(k)  (1) 

where  r(k)  is  a white  Gaussian  process  and  where  q(k)  is  a binary 
process,  which  can  only  take  on  the  values  0 or  1.  When  a reflection 
arrives  at  time  k we  set  q(k)«l,  and  when  no  reflection  arrives  at 
time  k we  set  q(k)"0.  Therefore  the  process  q(k)  is  determined  entirely 
by  the  locations  of  the  reflections,  so  that  the  magnitude,  information 
must  be  contained  in  the  process  r(k) . 

The  statistical  distribution  of  the  r (k)  is  completely  described  by 
its  variance,  C,  since  its  mean  is  zero.  The  distribution  of  the  q(k) 
is  given  by 
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when  qOO-l 
when  q(k)H) 
otherwise 


Pr{q(k)} 


where  the  paraaeter  A can  be  thought  of  as  the  average  ixafcer  of 
reflections  per  sample.  Using  the  product  model,  we  intend  to  find 
max  imun- likelihood  estimates  for  both  the  q(k)  and  the  r (k) , then  use 
the  invariance  property  of  maximun-likelihood  estimates  to  obtain 
estimates  of  the  ti(k)  using 


For  the  case  when  the  reflection  locations  are  known,  so  that  the 


using  an  optimal  another  by  incorporating  the  q(k)  into  the  state 
vector  model.  This  results  from  r(k)  being  a purely  Gaussian  process, 
and  the  maximum-likelihood  estimates  of  a Gaussian  process  are  the  same 
as  the  min  Jama-variance  estimates.  Figure  5 demonstrates  how  the  q(k) 
enter  into  the  convolution  model.  Figure  6 shows  the  resulting 
estimates  for  the  case  when  the  reflection  locations  are  known.  As  one 


can  see,  these  estimates  are  far  superior  to  those  obtained  without  such 
a priori  knowledge,  which  were  depicted  in  Figure  1. 

Another  way  to  estimate  the  v(k)  when  the  reflection  locations  are 
know*  is  to  incorporate  the  q(k)  into  the  variance  model.  That  is. 


E{y2(k) } - C 


which  was  used  to  obtain  the  minimun-variance  estimates,  we  could  use 
the  time-varying  conditional  variance 


E{y2(k)  | q(k) } » C q(k) 


and  model  the  y(k)  as  being  the  input  to  the  convolution  model.  Hiis 
way  the  optimal  smoother  will  estimate  y(k)  directly. 

A 

For  the  case  when  only  the  max imun-likel ihood  estimates  q(k)  are 
known  a priori,  the  maximum-likelihood  estimates  of  the  r(k)  can  be 

A 

obtained  using  an  optimal  anoother  by  incorporating  the  q(k)  into  the 
state  vector  model,  in  place  of  the  now  unknown  q(k) . This  can  be  shown 
by  factoring  the  likelihood  expression  into  a term  which  depends  on  the 
r(k)  and  a term  which  depends  of  the  q(k) . The  term  which  depends  of 
the  r(k)  will  have  the  form  of  a conditional  Gaussian  density  function, 

A 

but  where  q(k)  now  appears  in  place  of  q(k) . Therefore  the  problem  of 
estimating  the  reflection  amplitudes  has  a known  closed  form  solution, 
but  it  requires  first  estimating  the  reflection  locations. 


must  be  estimated  by  maximising  the  likelihood  expression  , which  is 
given  by 


Mg*’  | .«"»!  - pd'"1  | a(,W)  Pt(g,MI) 
using  the  notation 


(6) 


*(M)  - ts(l),  s(2),  ...,  x(M)]* 


(7) 


and 

(mi 

£ - lq(0),  q(l),  ...,  q<M-l)l'  . (8) 

Given  the  q(k) , we  can  express  the  x(k)  as  linear  combinations  of  the 
r(k)  plus  the  recording  noise,  both  of  which  are  modeled  as  Gaussian; 
hence,  the  probability  density  function  p{z(M)  |cj(M) } is  simply  a 
multivariate  Gaussian  density  function.  Since  the  reflection  locations 
are  statistically  independent,  the  probability  distribution  of  the  q(k) 
is  given  by 


- JJJ  Pr{q(k)}  (9) 

where  Pr{q(k) } was  given  by  equation  (2) . 

Na  can  compute  the  likelihood  for  any  given  set  of  observations 
tW  «d  estimates  £(M)  efficiently  using  a Kalman  filter  by  incor- 
porating the  estimates  q(k)  into  the  model,  since  the  likelihood  is  a 


simple  function  of  the  resulting  innovations  process  and  variance  [5] 

But  while  we  can  compute  the  likelihood  easily,  it  is  impossible  to 

(M) 

maximize  it  with  respect  to  the  unknown  vector  q in  any  iasonable 


Since  there  are  non-linear  interactions  between  the  q(k)  in  the 
likelihood  expression,  we  cannot  decouple  the  maximization  with  respect 
to  the  individual  q(k) , but  rather  must  maximize  with  respect  to  the 
entire  sequence.  Since  the  q(k)  are  discrete  valued,  maximization  will 
consist  of  comparing  the  likelihoods  of  different  sequences  (rather  than 
setting  a gradient  to  zero) . Finally,  since  each  q(k)  can  take  on  two 
possible  values  and  since  there  are  M samples  in  the  sequence,  there  is 

H (H) 

a total  of  2 possible  sequences  £ for  which  one  must  compare  the 
likelihoods.  When  M is  on  the  order  of  several  hundred,  this  would 


While  we  may  not  be  able  to  determine  the  global  maximum  of  the 
likelihood  expression,  we  can  always  design  a method  for  detecting 
significant  reflections;  that  is,  we  can  determine  relatively  high 
likelihood  or  near  maximum  likelihood  estimates  for  the  q(k) . Further 
since  we  can  easily  compute  the  likelihood  for  a given  estimated 

* (H) 

sequence  £ , we  can  compare  the  relative  performance  of  different 

detectors.  We  have  found  one  such  detector  which  performs  very  well 


This  detector  uses  an  approximation  to  the  likelihood  ratio  for 


q(k)  of  the  form 


A(k) 


l{  £(k+A)(k)  3 q(k)-l  | z(k+A)} 


L{  £(k+A)(k)  3 q(k)-0  I z(k+*)> 


(10) 


where  I»{  • | * } is  given  by  equation  (6)  and  where  the  k+t  length  vector 
a(k+t) (k)  is  defined  as 


fl(k+A) (k)  - [q(0) , q(l) , ...»  q(k-l),  q(k) , \,  ...,  x] ’ . (11) 


This  ratio  compares  the  likelihood  for  two  particular  sequences  which 
differ  only  in  their  values  for  q(k) . Further,  these  likelihoods  use 
only  k+i  observations,  where  £ is  some  small  integer,  rather  than  using 
all  M available  observations.  These  sequences  are  generated  by 

A A A 

recursively  detecting  q(0) , q(l),  ...  q(M)  using 


q(k) 


{; 


when  A(k)>l 
0 when  A(k)<l 


(12) 


A * fV+i.1 

for  k-0,l...M-l.  The  value  of  q(j)  in  £ (k)  for  j>k  is  replaced  by 


the  expected  value. 


E{q(j)}  - X 


(13) 


Figure  7 depicts  the  estimates  which  result  from  this  detector 


- 8 - 


using  £*5.  The  reflection  amplitudes  were  estimated  by  incorporating 

A 

tiie  detected  q(k)  into  the  state  vector  model.  These  estimates  are 
obviously  not  as  good  as  those  in  Figure  6,  in  spite  of  the  fact  that 
their  likelihood  is  many  times  higher.  The  reason  for  this  is  that  the 
seven  missed  reflections  were  not  large  enough  to  be  likely.  Further, 

A /u\ 

we  have  determined  that  this  estimated  sequence  q is  more  likely  than 
any  of  its  nearest  neighbors. 

Because  of  common  terms  in  the  two  sequences  used  in  equation  (10) , 
this  approximate  likelihood  ratio  can  be  computed  by  running  two  Kalman 
filters  for  only  i samples  each,  so  that  detecting  the  entire  sequence 

* /u\ 

3 requires  about  2i  times  as  much  computer  time  as  computing  its 
likelihood.  While  increasing  l increases  the  computational  burden,  it 
also  improves  the  resulting  likelihood.  But  since  the  marginal 
improvement  from  increasing  l drops  off  rapidly  as  l increases,  a 
relatively  snail  i will  generally  suffice.  Figure  8 demonstrates  the 
effect  of  varying  l . 

As  shown  in  Figure  9,  as  the  signal- to- noise  ratio  decreases,  the 
nunber  of  missed  detections  increases. 


VI.  Future  work 

If  one  is  more  interested  in  reducing  the  nunber  of  missed 
detections  than  in  maintaining  the  invariance  property,  then  one  might 
prefer  to  obtain  Bayesian  estimates  of  the  reflection  locations.  The 
solution  to  the  Bayesian  detection  problem  consists  of  comparing  the 


I 

1 


likelihood  ratio  to  some  threshold  other  than  unity  [6] . Since  our 
suboptimal  detector  computes  an  approximate  likelihood  ratio,  one  can 
easily  extend  this  approach  to  obtain  approximate  Bayesian  estimates. 

One  advantage  of  the  maximum- likelihood  approach  is  that  it  can 
handle  the  case  when  the  source  wavelet,  or  any  other  model  parameter, 
is  unknown.  This  is  done  by  obtaining  maximun-likel ihood  estimates  for 
these  unknown  parameters,  which  is  a conventional  non-linear  optimiz- 
ation problem.  We  have  shown  that  this  approach  can  resolve  the  phase 
ambiguity  problem  in  wavelet  extraction.  This  is  the  area  of  our 
present  interest. 

Also  this  approach  can  be  extended  to  include  any  bimodal  or 
multimodal  Gaussian  process.  For  example,  the  Earth  model  may  include 
snail  reflectors  within  the  layers,  so  that  the  reflection  coefficient 
sequence  is  actually  drawn  from  two  distributions:  one  for  the  layer 
interfaces  and  one  for  the  intra-layer  structures.  If  we  model  these 
magnitude  distributions  as  being  Gaussian  with  differing  variances,  then 
we  can  represent  the  reflection  coefficient  sequence  using  our  product 
model  by  letting  q(k)  take  on  two  non-zero  values. 
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VII.  Conclusion 

~~ 

By  representing  the  reflection  coefficient  sequence  as  a product  of 
a Gaussian  process  and  a binary  process,  the  estimation  problem  is 
separated  into  two  simpler  problems.  First,  the  significant  reflections 
are  detected.  Second,  the  amplitudes  are  estimated  using  an  optimal 
smoother  which  incorporates  the  solution  to  the  first  problem  in  the 
state  vector  model. 

While  it  is  impossible  to  definitely  determine  the  reflection 
locations  which  maximize  the  likelihood  in  a reasonable  amount  of  time, 
we  have  presented  an  efficient  method  for  recursively  detecting  the 
significant  reflections  using  an  approximation  to  the  likelihood  ratio. 
This  subopt imal  detector  has  performed  very  well  on  synthetic  data. 
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CONVOLUTION  MODEL 


>z(k) 


fiW 


/*00  is  the  desired  reflection  sequence. 
z(k)  is  the  observed  data. 


Pig.  2.  Convolution  Model 
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ESTIMATING  AMPLITUDES  GIVEN  LOCATIONS 

When  the  q(k)  ore  known,  the  maximum* likelihood 
estimates  for  the  r(k)  can  be  obtained  by  using  on 
optimal  smoother  where  the  q(k)  are  included  in  the 
state  vector  model  : 


r(k) 


We  con  then  estimate  fi[k)  using 

ft(k)  * *(k)  q (k) 


Fig.  5.  Incorporating  q(k)  into  tha  Convolution  Modal 
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Fig.  7.  Detected  Estimates,  SNR-10 
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Fig.  8.  Log  Likelihood  vs.  I 
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